Local Metabolism Plots

# set the root to the project directory, not the rmd directory
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# default to not printing the code
knitr::opts_chunk$set(echo = FALSE)

# libraries
library(tidyverse)

package 㤼㸱tidyverse㤼㸲 was built under R version 4.0.5Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages ------------------------------------------------------------ tidyverse 1.3.1 --
v ggplot2 3.3.3     v purrr   0.3.4
v tibble  3.1.2     v dplyr   1.0.6
v tidyr   1.1.3     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.1
package 㤼㸱ggplot2㤼㸲 was built under R version 4.0.3package 㤼㸱tibble㤼㸲 was built under R version 4.0.5package 㤼㸱tidyr㤼㸲 was built under R version 4.0.5package 㤼㸱readr㤼㸲 was built under R version 4.0.3package 㤼㸱dplyr㤼㸲 was built under R version 4.0.5package 㤼㸱forcats㤼㸲 was built under R version 4.0.3-- Conflicts --------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(sf)

Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

library(stars)

package 㤼㸱stars㤼㸲 was built under R version 4.0.5Loading required package: abind

library(tmap)

package 㤼㸱tmap㤼㸲 was built under R version 4.0.5Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(transformr) # not needed for this doc at present?

package 㤼㸱transformr㤼㸲 was built under R version 4.0.5
Attaching package: 㤼㸱transformr㤼㸲

The following object is masked from 㤼㸱package:sf㤼㸲:

    st_normalize

library(gganimate) # not needed here at present

package 㤼㸱gganimate㤼㸲 was built under R version 4.0.5

library(viridis)

package 㤼㸱viridis㤼㸲 was built under R version 4.0.5Loading required package: viridisLite
package 㤼㸱viridisLite㤼㸲 was built under R version 4.0.5

library(colorspace)

package 㤼㸱colorspace㤼㸲 was built under R version 4.0.5

# source in the data management

# source('directorySet.R')
source(file.path('Scripts', 'Scenarios', 'plotting', 'metabolismLocalAndStatic.R'))

          sysname           release           version          nodename           machine 
        "Windows"          "10 x64"     "build 19042" "DESKTOP-NU5UDHD"          "x86-64" 
            login              user    effective_user 
          "Galen"           "Galen"           "Galen" 

Purpose here is primarily to go through some plots, look at them, and make static versions for the paper. This is mostly a copy-paste and translate over from metabolismLocalAndStatic. I’ve done this to split out the data read-in, function definitions, and plotting (and particularly the plot trying-out). I’ve been using these notebooks for plot testing, because it’s really nice to be able to scroll around and see what we’ve done, rather than have to keep re-making plot objects when we forget what they look like.

The catch with this approach is the quick and dirty prototyping isn’t quite so quick, and this document ends up taking FOREVER to load. Might switch back? But it sure is nice to actually see what each figure I made looks like, especially when we’re likely to drop this for extended periods.

Set the date for examples

availDays <- st_get_dimension_values(weraiCropTemp, which = 'time')
datewanted <- as.character(availDays[17])

Static versions with drivers and outcomes

Make a static version using the tmap functions

suppressMessages(tmap_mode('plot'))

allfun(tempObj = weraiCropTemp, tempAtt = 1,
                   inunObj = weraiCropInun, inunAtt = 1,
                   gppObj = weraiCropPredGPP, gppAtt = 1,
                   erObj = weraiCropPredER, erAtt = 1,
                   datewanted)

Clean that up just a bit, with a grey background. probably pull the title off but it’s nice to have as reference while i’m making these

tm_grid_static <- tmap_arrange(tempfun(weraiCropTemp, 1, datewanted) +
                                 tm_layout(bg.color = "grey85"),
                               inunfun(weraiCropInun, 1, datewanted, titled = FALSE) +
                                 tm_layout(bg.color = "grey85"),
                               gppfun(weraiCropPredGPP, 1, datewanted, titled = FALSE) +
                                 tm_layout(bg.color = "grey85"),
                               erfun(weraiCropPredER, 1, datewanted, titled = FALSE) +
                                 tm_layout(bg.color = "grey85"))
tm_grid_static

I think for completeness let’s make a ggplot version too and see which is better. Takes more work, but have more control (probably just because I speak ggplot)

gg_grid_static <- ggpubr::ggarrange(
  tempfun(weraiCropTemp, 1, datewanted, plotPkg = 'ggplot') + 
    guides(fill = guide_legend(title.position = 'top')) +
    theme_grey(base_size = 8) + 
    theme(legend.position = 'bottom', 
          legend.background = element_blank(),
          legend.key.size = unit(0.3, 'cm')),
  inunfun(weraiCropInun, 1, datewanted, titled = FALSE, plotPkg = 'ggplot') + 
    guides(fill = guide_legend(title.position = 'top')) +
    theme_grey(base_size = 8) + 
    theme(legend.position = 'bottom', 
          legend.background = element_blank(),
          legend.key.size = unit(0.3, 'cm')),
  gppfun(weraiCropPredGPP, 1, datewanted, titled = FALSE, plotPkg = 'ggplot') + 
    guides(fill = guide_legend(title.position = 'top')) +
    theme_grey(base_size = 8) + 
    theme(legend.position = 'bottom', 
          legend.background = element_blank(),
          legend.key.size = unit(0.3, 'cm')),
  erfun(weraiCropPredER, 1, datewanted, titled = FALSE, plotPkg = 'ggplot') + 
    guides(fill = guide_legend(title.position = 'top')) +
    theme_grey(base_size = 8) + 
    theme(legend.position = 'bottom', 
          legend.background = element_blank(),
          legend.key.size = unit(0.3, 'cm')),
  ncol = 2, nrow = 2)
print(gg_grid_static)

Code block to print and make those into figs

pdf(file.path(scriptOut, 'Werai_tm.pdf'), 
    onefile = FALSE, height = 12/2.54, width = 16/2.54, useDingbats = FALSE)
print(tm_grid_static)
dev.off()

null device 
          1 

png(file.path(scriptOut, 'Werai_tm.png'), 
    height = 12/2.54, width = 16/2.54, units = 'in', res = 300)

print(tm_grid_static)
dev.off()

null device 
          1 

# Can I print those?
pdf(file.path(scriptOut, 'Werai_gg.pdf'), 
    onefile = FALSE, height = 12/2.54, width = 16/2.54, useDingbats = FALSE)

print(gg_grid_static)
dev.off()

null device 
          1 

png(file.path(scriptOut, 'Werai_gg.png'), 
    height = 12/2.54, width = 16/2.54, units = 'in', res = 300)

print(gg_grid_static)
dev.off()

null device 
          1 

Uncertainty

I think now that I have the functions, there’s a slicker way to do this, but I don’t want to reinvent the wheel at this point. Might want to turn this into a function?

ER

# ER
er_sf <- weraiCropPredER[1,,] %>% 
  st_as_sf() %>% 
  select(all_of(datewanted)) %>%
  rename(ER = 1) %>%
  mutate(logER = log10(1+ER), type = 'estimate')
er_sfPU <- weraiCropPredER[2,,] %>% 
  st_as_sf() %>% 
  select(all_of(datewanted)) %>%
  rename(ER = 1) %>%
  mutate(logER = log10(1+ER), type = 'upper')
er_sfPL <- weraiCropPredER[3,,] %>% 
  st_as_sf() %>% 
  select(all_of(datewanted)) %>%
  rename(ER = 1) %>%
  mutate(logER = log10(1+ER), type = 'lower')

er_sf_LMU <- bind_rows(er_sf, er_sfPU, er_sfPL) %>%
  mutate(namedPI = ifelse(type == 'lower', "Lower 95% PI",
                          ifelse(type == 'estimate', "Estimate",
                                 ifelse(type == 'upper', "Upper 95% PI", 
                                        'SCREWUP'))))

erlabel <- 'ER (kg 02/day)\nat max extent'

# Can get away with this because pretty ranges work out OK
erControl <- ersetup(data = weraiCropPredER, attnum = 1)

  er_gg_uncertain <- ggplot() +
    geom_sf(data = er_sf_LMU, mapping = aes(fill = logER), color = NA) +
    # geom_sf_label(data = ltimNoNorth, mapping = aes(label = ValleyName)) +
    coord_sf() +
    # Closest to the tmap
    scale_fill_stepsn(colors = erControl$erpal,
                      breaks = erControl$erbreaks[2:length(erControl$erbreaks)],
                      limits = c(min(erControl$erbreaks), max(erControl$erbreaks)),
                      labels = erControl$erlabels,
                      guide = 'legend',
                      name = erlabel) +
    facet_grid(fct_reorder(namedPI, .x = logER, .fun = max)~.) # lower, estimat and upper should always fall in that order for their maxes

er_gg_uncertain

GPP

# GPP
# I'm sure thgppe's a slick way to select all three attributes and do this in one
# go, but I just need to get this done
gpp_sf <- weraiCropPredGPP[1,,] %>% 
  st_as_sf() %>% 
  select(all_of(datewanted)) %>%
  rename(GPP = 1) %>%
  mutate(logGPP = log10(1+GPP), type = 'estimate')
gpp_sfPU <- weraiCropPredGPP[2,,] %>% 
  st_as_sf() %>% 
  select(all_of(datewanted)) %>%
  rename(GPP = 1) %>%
  mutate(logGPP = log10(1+GPP), type = 'upper')
gpp_sfPL <- weraiCropPredGPP[3,,] %>% 
  st_as_sf() %>% 
  select(all_of(datewanted)) %>%
  rename(GPP = 1) %>%
  mutate(logGPP = log10(1+GPP), type = 'lower')

gpp_sf_LMU <- bind_rows(gpp_sf, gpp_sfPU, gpp_sfPL) %>%
  mutate(namedPI = ifelse(type == 'lower', "Lower 95% PI",
                          ifelse(type == 'estimate', "Estimate",
                                 ifelse(type == 'upper', "Upper 95% PI", 
                                        'SCREWUP'))))

gpplabel <- 'GPP (kg 02/day)\nat max extent'

gppControl <- gppsetup(weraiCropPredGPP, attnum  = 1)

gpp_gg_uncertain <- ggplot() +
  geom_sf(data = gpp_sf_LMU, mapping = aes(fill = logGPP), color = NA) +
  # geom_sf_label(data = ltimNoNorth, mapping = aes(label = ValleyName)) +
  coord_sf() +
  # Closest to the tmap
  scale_fill_stepsn(colors = gppControl$gpppal,
                    breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                    limits = c(min(gppControl$gppbreaks), max(gppControl$gppbreaks)),
                    labels = gppControl$gpplabels,
                    guide = 'legend',
                    name = gpplabel) +
  facet_grid(fct_reorder(namedPI, .x = logGPP, .fun = max)~.) # lower, estimat and uppgpp should always fall in that ordgpp for their maxes

gpp_gg_uncertain

Combine the ER and GPP and save as figures

both_uncertain <- ggpubr::ggarrange(gpp_gg_uncertain + 
                                   guides(fill = guide_legend(title.position = 'top')) +
                                   theme_grey(base_size = 8) + 
                                   theme(legend.position = 'bottom', 
                                         legend.background = element_blank(),
                                         legend.key.size = unit(0.3, 'cm')), 
                                 er_gg_uncertain + 
                                   guides(fill = guide_legend(title.position = 'top')) +
                                   theme_grey(base_size = 8) + 
                                   theme(legend.position = 'bottom', 
                                         legend.background = element_blank(),
                                         legend.key.size = unit(0.3, 'cm')),
                                 ncol = 2, nrow = 1)
both_uncertain

# Just print the ggplots
# Can I print those?
pdf(file.path(scriptOut, 'Werai_uncertainty.pdf'), 
    onefile = FALSE, height = 12/2.54, width = 16/2.54, useDingbats = FALSE)

print(both_uncertain)
dev.off()

png 
  2 

png(file.path(scriptOut, 'Werai_uncertainty.png'), 
    height = 12/2.54, width = 16/2.54, units = 'in', res = 300)
print(both_uncertain)

dev.off()

png 
  2 

Bar charts

make a wide version of the data so I can have mins and maxes OR, should I do all of them and fct_order them?

Joining seems safe, but geographic joins are a mess. It doesn’t work at all. binding cols is actually safer and works

GPP

gpp_sf_LMU_wide <- bind_cols(gpp_sf[,-4], 
                           st_drop_geometry(rename(gpp_sfPU[,-c(4)], 
                                  UGPP = GPP, 
                                  upper_logGPP = logGPP)),
                           st_drop_geometry(rename(gpp_sfPL[,-4], 
                                  LGPP = GPP,
                                  lower_logGPP = logGPP))) %>%
  mutate(wetlandID = row_number())

Maybe for a subset of wetlands?

# Grab some set of wetlands
which(gpp_sf$geometry == gpp_sf_LMU$geometry[3])

[1] 3

exampleWets <- c(1:100)

bargpp <- ggplot(gpp_sf_LMU_wide[exampleWets, ], 
                 aes(x = wetlandID, y = logGPP, fill = logGPP)) +
  geom_col() + # the geom_bar equivalent when it's not counting cases.
  geom_errorbar(mapping = aes(ymin = lower_logGPP, ymax = upper_logGPP)) +
  scale_fill_stepsn(colors = gppControl$gpppal,
                    breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                    limits = c(min(gppControl$gppbreaks), max(gppControl$gppbreaks)),
                    labels = gppControl$gpplabels,
                    guide = 'legend',
                    name = gpplabel) +
  scale_y_continuous(breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                     labels = round(10^(gppControl$gppbreaks[2:length(gppControl$gppbreaks)]))) +
  labs(x = 'Wetland', y = gpplabel)
bargpp

I assume it is a mes to just throw them all on actually, with a fair amount of monkeying with the look, that’s OK

bargppall <- ggplot(gpp_sf_LMU_wide, 
                 aes(x = fct_reorder(as.factor(wetlandID), logGPP, .desc = TRUE), y = logGPP, fill = logGPP)) +
  geom_col() + # the geom_bar equivalent when it's not counting cases.
  geom_linerange(mapping = aes(ymin = lower_logGPP, ymax = upper_logGPP),
                 color = 'grey50', alpha = 0.2) +
  scale_fill_stepsn(colors = gppControl$gpppal,
                    breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                    limits = c(min(gppControl$gppbreaks), max(gppControl$gppbreaks)),
                    labels = gppControl$gpplabels,
                    guide = 'legend',
                    name = gpplabel) +
  scale_y_continuous(breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                     labels = round(10^(gppControl$gppbreaks[2:length(gppControl$gppbreaks)]))) +
  labs(x = 'Wetland', y = gpplabel)
bargppall

I actually like that that encompasses an area that integrates to total production across the Werai

I can do points as below, but the integration implied above is nice

bargppallpoint <- ggplot(gpp_sf_LMU_wide, 
                         aes(x = fct_reorder(as.factor(wetlandID), logGPP), 
                             y = logGPP, color = logGPP)) +
  geom_linerange(mapping = aes(ymin = lower_logGPP, ymax = upper_logGPP), color = 'grey50') +
  geom_point() + # the geom_bar equivalent when it's not counting cases.
  scale_color_stepsn(colors = gppControl$gpppal,
                     breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                     limits = c(min(gppControl$gppbreaks), max(gppControl$gppbreaks)),
                     labels = gppControl$gpplabels,
                     guide = 'legend',
                     name = gpplabel) +
  scale_y_continuous(breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                     labels = round(10^(gppControl$gppbreaks[2:length(gppControl$gppbreaks)]))) +
  labs(x = 'Wetland', y = gpplabel)
bargppallpoint

ER

Setup block. Again, this could be sorted with a function, but right now just copy-pasting

erlabel <- 'GPP (kg 02/day)\nat max extent'

er_sf_LMU_wide <- bind_cols(er_sf[,-4], 
                             st_drop_geometry(rename(er_sfPU[,-c(4)], 
                                                     UER = ER, 
                                                     upper_logER = logER)),
                             st_drop_geometry(rename(er_sfPL[,-4], 
                                                     LER = ER,
                                                     lower_logER = logER))) %>%
  mutate(wetlandID = row_number())

Barplot

barerall <- ggplot(er_sf_LMU_wide, 
                    aes(x = fct_reorder(as.factor(wetlandID), logER, .desc = TRUE), y = logER, fill = logER)) +
  geom_col() + # the geom_bar equivalent when it's not counting cases.
  geom_linerange(mapping = aes(ymin = lower_logER, ymax = upper_logER),
                 color = 'grey50', alpha = 0.2) +
  scale_fill_stepsn(colors = erControl$erpal,
                    breaks = erControl$erbreaks[2:length(erControl$erbreaks)],
                    limits = c(min(erControl$erbreaks), max(erControl$erbreaks)),
                    labels = erControl$erlabels,
                    guide = 'legend',
                    name = erlabel) +
  scale_y_continuous(breaks = erControl$erbreaks[2:length(erControl$erbreaks)],
                     labels = round(10^(erControl$erbreaks[2:length(erControl$erbreaks)]))) +
  labs(x = 'Wetland', y = erlabel)
barerall

Both together, with GPP up and ER down

Setup

# Can I make one with them both there going opposite direction?
both_sf_LMU_wide <- bind_cols(mutate(gpp_sf_LMU_wide), 
                              mutate(st_drop_geometry(er_sf_LMU_wide[,-c(8)])))

# Need new labels and breaks. Could cobble, but getting confusing and misaligned
# bothbreaks_log <- c(-rev(erbreaks_log[2:length(erbreaks_log)]), gppbreaks_log)
 
# # and those breaks might not quite yield 10, so maximise the palette differences
# erpal_log <- sequential_hcl(length(erbreaks_log)-1, palette = 'Purples', rev = TRUE)
# I think here I want continuous colors. but maybe still a broken up legend?

# have to do the negatives before delogging
# bothlabels_log <- c(-10^rev(erbreaks_log), 10^gppbreaks_log)
# 
# bothlabels_log <- format(bothlabels_log, big.mark=",", 
                       # scientific=FALSE, trim = TRUE, digits = 0)

# This is really getting complex. Why don't I just use 0,1,10,100,1000,10000?
# They've both been shifted by 1 to log, so
# This is just used to get the labels, the actual log-scale is plotted and so the zeros are done correctly
bothbreaks_log <- c(-10^(4:0), 10^(1:4))
bothbreaks_log[5] <- 0 # Because both were shifted so 1 = 0
bothlabels_log <- as.character(bothbreaks_log)


# erstart <- erlabels_log[1:(length(erlabels_log)-1)]
# erstart[1] <- "0" # instead of 1
# bothlabels_log <- paste0(erstart, ' to ', erlabels_log[2:length(erlabels_log)])
# erlabels_log

# AAAA try to get some lables for the x. I JUST WANT IT TO HAVE THE NUMBER BUT IT WON"T DO IT
reordx <- fct_reorder(as.factor(both_sf_LMU_wide$wetlandID), 
                      both_sf_LMU_wide$logGPP, .desc = TRUE)
# I'm confused. that didn't work. this is anoyuting

Try a plot

barboth <- ggplot(both_sf_LMU_wide) +
  geom_col(mapping = aes(x = fct_reorder(as.factor(wetlandID), 
                                         logER, .desc = TRUE), 
                         y = -logER, fill = -logER)) + # the geom_bar equivalent when it's not counting cases.
  geom_linerange(mapping = aes(x = fct_reorder(as.factor(wetlandID), 
                                               logER, .desc = TRUE),
                               ymin = -lower_logER, ymax = -upper_logER),
                 color = 'grey50', alpha = 0.2) +
  # scale_fill_stepsn(colors = erpal_log,
  #                   breaks = erbreaks_log[2:length(erbreaks_log)],
  #                   limits = c(min(erbreaks_log), max(erbreaks_log)),
  #                   labels = erlabels_log,
  #                   guide = 'legend',
  #                   name = erlabel) +
 
  # scale_fill_stepsn(colors = c(rev(erpal_log), gpppal_log),
  #            breaks = c(-rev(erbreaks_log[2:length(erbreaks_log)]), 0,
  #                       gppbreaks_log[2:length(gppbreaks_log)]),
  #            limits = c(min(-erbreaks_log), max(gppbreaks_log)),
  #            # labels = c(erlabels_log, gpplabels_log),
  #            guide = 'legend',
  #            name = 'kg O2\nat max extent') +
  geom_col(mapping = aes(x = fct_reorder(as.factor(wetlandID), 
                                         logGPP, .desc = TRUE), 
                         y = logGPP, fill = logGPP)) +
  geom_linerange(mapping = aes(x = fct_reorder(as.factor(wetlandID), 
                                               logGPP, .desc = TRUE),
                               ymin = lower_logGPP, ymax = upper_logGPP),
                 color = 'grey50', alpha = 0.2) +
  geom_line(mapping = aes(x = fct_reorder(as.factor(wetlandID), 
                                           logGPP, .desc = TRUE), 
                           y = logGPP-logER, group = NA)) +
  
  # This palette is very slightly different because it's not binned and GPP above uses Emerald instead of green, but it's sure close and way easier
  scale_fill_continuous_diverging(palette = "Purple-Green",
                                  limits = c(-4, 4), 
                                  breaks = -4:4,
                                  labels = bothlabels_log) +
  scale_y_continuous(limits = c(-4, 4),
                     breaks = -4:4,
                     labels = bothlabels_log) +
  # scale_y_continuous(breaks = c(-erbreaks_log[2:length(erbreaks_log)], 
  #                               gppbreaks_log[2:length(gppbreaks_log)]),
  #                    labels = c(-round(10^(erbreaks_log[2:length(erbreaks_log)])), 
  #                               round(10^(gppbreaks_log[2:length(gppbreaks_log)])))) +
  labs(x = 'Wetland (ordered by GPP)', y = 'kg O2', fill = 'kg O2') +
  # argh doesn't work because has been reordered
  # scale_x_discrete(breaks = as.character(c(1, 250, 500, 750, 1000))) +
  theme(legend.position = c(0.75, 0.54)) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = c(0.75, 0.55))
barboth

The inability to label x in a reasonable way is infuriating. Can I do it differently? Try making a labelled dataset that’s arranged how I was doing fct_reorder.

both_sf_LMU_wide2 <- both_sf_LMU_wide %>% 
  arrange(desc(logGPP)) %>% 
  mutate(GPPrank = row_number())

Make the plot. Interesting that the notebook plotter does a worse job than the plots panel for a script. The funny gaps aren’t really there- to see, run barboth in the console.

barboth <- ggplot(both_sf_LMU_wide2) +
  geom_col(mapping = aes(x = GPPrank, 
                         y = -logER, fill = -logER)) + # the geom_bar equivalent when it's not counting cases.
  geom_linerange(mapping = aes(x = GPPrank,
                               ymin = -lower_logER, ymax = -upper_logER),
                 color = 'grey50', alpha = 0.2) +

geom_col(mapping = aes(x = GPPrank, 
                       y = logGPP, fill = logGPP)) +
  geom_linerange(mapping = aes(x = GPPrank,
                               ymin = lower_logGPP, ymax = upper_logGPP),
                 color = 'grey50', alpha = 0.2) +
  geom_line(mapping = aes(x = GPPrank, 
                          y = logGPP-logER, group = NA)) +
  
  # This palette is very slightly different because it's not binned and GPP above uses Emerald instead of green, but it's sure close and way easier
  scale_fill_continuous_diverging(palette = "Purple-Green",
                                  limits = c(-4, 4), 
                                  breaks = -4:4,
                                  labels = bothlabels_log) +
  scale_y_continuous(limits = c(-4, 4),
                     breaks = -4:4,
                     labels = bothlabels_log) +
  labs(x = 'Wetland (ordered by GPP)', y = 'kg O2', fill = 'kg O2') +
  # scale_x_discrete(breaks = c(1, 250, 500, 750, 1000)) +
  theme_bw(base_size = 8) + theme(legend.position = c(0.75, 0.53),
                                  panel.grid.major = element_blank(),
                                  panel.grid.minor = element_blank())
barboth

Save that as a figure

pdf(file.path(scriptOut, 'Werai_mirrorgram.pdf'), 
    onefile = FALSE, height = 8/2.54, width = 12/2.54, useDingbats = FALSE)
print(barboth)
dev.off()

null device 
          1 

png(file.path(scriptOut, 'Werai_mirrorgram.png'), 
    height = 8/2.54, width = 12/2.54, units = 'in', res = 300)

print(barboth)
dev.off()

null device 
          1 

What about Darren’s fingerprints? They would maybe take the GPP and ER and make a density with them? OUt of what? The wetlands? Could work. Would need to sort out the weighting by volume. it’s already in there, but that means it will give each location the same weight. I think a better thing for the fingerprints would be to remove it, and then give each wetland its predicted per liter rate and weight by volume IE use the straight predictions and then weight by inundaton Should be straightforward, but likely won’t be

Would be a really cool way to look at scenarios. especially if I can do it for whole catchments

BUT, because GPP and ER are both linear fits of temp, they will exactly predict each other according to a linear relationship, and so there won’t be any 2-d variation and all the points will fall on a line. IE for a given GPP there will only be one ER, and so the fingerprint will be a line. Cool idea, but we’d need more info about their variance relative to each other

Annual reporting

What about making some examples to illustrate annual (or arbitrary time period) reporting?

First, establish the time periods and do the aggregations and other data organisation

interDates <- as.POSIXct(c("2014-06-30", "2015-06-30", "2016-06-30", "2017-06-30", "2018-06-30", "2019-06-30"))

tempannual <- tempaggregate(weraiCropTemp, by_t = as.Date(interDates), FUN = mean, na.rm = TRUE)
inunannual <- tempaggregate(weraiCropInun, by_t = as.Date(interDates), FUN = sum, na.rm = TRUE)
gppannual <- tempaggregate(weraiCropPredGPP[1,,], by_t = as.Date(interDates), FUN = sum, na.rm = TRUE)
erannual <- tempaggregate(weraiCropPredER[1,,], by_t = as.Date(interDates), FUN = sum, na.rm = TRUE)

# Data organisation
# AHHH the flopped dims
if (attributes(st_dimensions(gppannual))$name[1] != 'geometry') {
  gppannual <- aperm(gppannual, c(2,1))
}
if (attributes(st_dimensions(erannual))$name[1] != 'geometry') {
  erannual <- aperm(erannual, c(2,1))
}

# GPP
gppYear_sf <- gppannual %>% 
  st_as_sf() %>% 
  pivot_longer(cols = -geometry, names_to = 'WaterYear', values_to = 'GPP') %>%
  mutate(logGPP =  log10(1+GPP)) %>%
  st_as_sf() # REALLY???

number of items to replace is not a multiple of replacement length

gppYear_sf

Simple feature collection with 5815 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1086148 ymin: -3924954 xmax: 1148710 ymax: -3893940
Projected CRS: GDA94 / Australian Albers

max(gppYear_sf$logGPP) # is within 

[1] 3.412777

gppControl$gppbreaks

[1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

# ER
erYear_sf <- erannual %>% 
  st_as_sf() %>% 
  pivot_longer(cols = -geometry, names_to = 'WaterYear', values_to = 'ER') %>%
  mutate(logER =  log10(1+ER)) %>%
  st_as_sf() # REALLY???

number of items to replace is not a multiple of replacement length

erYear_sf

Simple feature collection with 5815 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1086148 ymin: -3924954 xmax: 1148710 ymax: -3893940
Projected CRS: GDA94 / Australian Albers

# the range should still work for tthe colors
max(erYear_sf$logER) # is within 

[1] 3.779558

erControl$erbreaks

[1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

I don’t think I use this, but I did make a full-werai aggregation

# aggregate to werai
# temp area-weighted
areas <- tempannual %>%
  st_geometry() %>%
  st_area() %>%
  as.numeric()
tempW <- catchAgg <- catchAggW(strict = tempannual, strictWeights = areas,
                               FUN = mean, summaryPoly = ramsarW1)
  
inunW <- aggregate(inunannual, 
                      by = ramsarW1, 
                      FUN = sum, na.rm = TRUE)

gppW <- aggregate(gppannual, 
                   by = ramsarW1, 
                   FUN = sum, na.rm = TRUE)
erW <- aggregate(erannual, 
                   by = ramsarW1, 
                   FUN = sum, na.rm = TRUE)

Tmap Plots

Make a facetted tmap

Again, can I do this by modifying what I pass to the plotting functions? Almost certainly

The units are weird, because this is added up across bimonths. It’s neither the total or the average.Not really sure what to call it

Change the tmap_mode() depending on the purpose. DO NOT CHANGE TO VIEW AND TRY TO RUN HERE. IF YOU WANT TO VIEW, CHANGE AND RUN IN CONSOLE. The markdown will run it, but incredibly slowly.

GPP

tmap_mode('plot')

tmap mode set to plotting

gppyrlabel <- 'Total Yearly GPP (kg 02)\nat max bimonth extents'

  gppAnnual_tm <- gppYear_sf %>%
    filter(WaterYear != '2014-06-29') %>% # because a 5-panel is ugly
    tm_shape() +
    tm_fill(col = 'logGPP',
            palette = gppControl$gpppal,
            breaks = gppControl$gppbreaks,
            labels = gppControl$gpplabels,
            title = gppyrlabel) + 
    tm_facets(by = 'WaterYear')
  gppAnnual_tm

ER

eryrlabel <- 'Total Yearly ER (kg 02)\nat max bimonth extents'
  
  erAnnual_tm <- erYear_sf %>%
    filter(WaterYear != '2014-06-29') %>% # because a 5-panel is ugly
    tm_shape() +
    tm_fill(col = 'logER',
            palette = erControl$erpal,
            breaks = erControl$erbreaks,
            labels = erControl$erlabels,
            title = eryrlabel) + 
    tm_facets(by = 'WaterYear')
  erAnnual_tm

GGPLOT

gpp

gppAnnual_gg <- gppYear_sf %>%
    filter(WaterYear != '2014-06-29') %>% # because a 5-panel is ugly
    ggplot() +
    geom_sf(mapping = aes(fill = logGPP), color = NA) +
    coord_sf() +
    # Closest to the tmap
    scale_fill_stepsn(colors = gppControl$gpppal,
                      breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                      limits = c(min(gppControl$gppbreaks), max(gppControl$gppbreaks)),
                      labels = gppControl$gpplabels,
                      guide = 'legend',
                      name = gppyrlabel) +
    facet_wrap(vars(WaterYear))
  gppAnnual_gg

ER

erAnnual_gg <- erYear_sf %>%
    filter(WaterYear != '2014-06-29') %>% # because a 5-panel is ugly
    ggplot() +
    geom_sf(mapping = aes(fill = logER), color = NA) +
    coord_sf() +
    # Closest to the tmap
    scale_fill_stepsn(colors = erControl$erpal,
                      breaks = erControl$erbreaks[2:length(erControl$erbreaks)],
                      limits = c(min(erControl$erbreaks), max(erControl$erbreaks)),
                      labels = erControl$erlabels,
                      guide = 'legend',
                      name = eryrlabel) +
    facet_wrap(vars(WaterYear))
  erAnnual_gg

Stack 3 years of gpp and ER together for the document

gppAnnual_gg3 <- gppYear_sf %>%
    filter(WaterYear != '2014-06-29' & WaterYear != '2018-06-29') %>% # because a 5-panel is ugly
    ggplot() +
    geom_sf(mapping = aes(fill = logGPP), color = NA) +
    coord_sf() +
    # Closest to the tmap
    scale_fill_stepsn(colors = gppControl$gpppal,
                      breaks = gppControl$gppbreaks[2:length(gppControl$gppbreaks)],
                      limits = c(min(gppControl$gppbreaks), max(gppControl$gppbreaks)),
                      labels = gppControl$gpplabels,
                      guide = 'legend',
                      name = gppyrlabel) +
    facet_grid(WaterYear~.) + 
    theme(legend.position = 'bottom')
  # gppAnnual_gg3 
  
  erAnnual_gg3 <- erYear_sf %>%
    filter(WaterYear != '2014-06-29' & WaterYear != '2018-06-29') %>% # because a 5-panel is ugly
    ggplot() +
    geom_sf(mapping = aes(fill = logER), color = NA) +
    coord_sf() +
    # Closest to the tmap
    scale_fill_stepsn(colors = erControl$erpal,
                      breaks = erControl$erbreaks[2:length(erControl$erbreaks)],
                      limits = c(min(erControl$erbreaks), max(erControl$erbreaks)),
                      labels = erControl$erlabels,
                      guide = 'legend',
                      name = eryrlabel) +
    facet_grid(WaterYear~.) + 
    theme(legend.position = 'bottom')
  # erAnnual_gg3

  annualGPPER <- ggpubr::ggarrange(gppAnnual_gg3 + 
                                     guides(fill = guide_legend(title.position = 'top')) +
                                     theme_grey(base_size = 8) + 
                                     theme(legend.position = 'bottom', 
                                           legend.background = element_blank(),
                                           legend.key.size = unit(0.3, 'cm')), 
                    erAnnual_gg3 + 
                      guides(fill = guide_legend(title.position = 'top')) +
                      theme_grey(base_size = 8) + 
                      theme(legend.position = 'bottom', 
                            legend.background = element_blank(),
                            legend.key.size = unit(0.3, 'cm')),
                    ncol = 2, nrow = 1)
  annualGPPER

Print out the ggplots for the doc

 # Just print the ggplots
  # Can I print those?
  pdf(file.path(scriptOut, 'Werai_gg_Annual.pdf'), 
      onefile = FALSE, height = 12/2.54, width = 16/2.54, useDingbats = FALSE)
  print(annualGPPER)
  dev.off()

null device 
          1 

  
  png(file.path(scriptOut, 'Werai_gg_Annual.png'), 
      height = 12/2.54, width = 16/2.54, units = 'in', res = 300)

  print(annualGPPER)
  dev.off()

null device 
          1